Accounting for uncertainty in simulation greenhouse gas emissions from dairy cattle in Kenya
Eike Luedeling, Cory Whitney, Todd Rosenstock
Set up
Install and load required packages
devtools::install_github("CWWhitney/uncertainty")
library(DiagrammeR)
library(decisionSupport)
library(ggplot2)
library(kableExtra)
library(magrittr)
library(plyr)
library(rmarkdown)
library(uncertainty)
library(chillR)
library(raster)Load emission calculation function
The formulas come from Chapter 10 of the 2006 IPCC Guidelines for National Greenhouse Gas Inventories (Dong et al. 2006) (Volume 4 - Agriculture, Forestry and Other Land Use). The document can be found on the IPCC website. The chapter is part of the book Good Practice Guidance and Uncertainty Management in National Greenhouse Gas Inventories (IPCC 2021), also on the the IPCC website. Many variables are based on the work done in Variation in the carbon footprint of milk production on smallholder dairy farms in central Kenya (Wilkes et al. 2020) and in Methods and guidance to support MRV of livestock emissions (Wilkes et al. 2019). We run the model with the decisionSupport package (Luedeling et al. 2022).
# Generate the GHG emissions function
source("GHG_emissions_function.R", local = knitr::knit_global())This equation takes a wide range of inputs. Many of these are usually not known precisely, yet the classic IPCC equation does not account for potential errors and uncertainties. Here, we express this uncertainty by defining distributions for each variable. The following table shows all the variables, as well as examples of distributions that may describe their quantities. Distributions are defined according to the distribution shape, as well as upper and lower bounds of 90% confidence intervals.
### # This project will be stored in relative file paths (i.e. same or lower folders)
input_table <- "data/livestock_ghg_input_table.csv"
options(knitr.kable.NA = '', encoding = "Windows-1252")
read.csv(input_table)[,c(1:2,4:6)] %>%
kableExtra::kbl(digits = 2) %>%
kableExtra::kable_classic_2(full_width = F)| Description | lower | upper | distribution | variable |
|---|---|---|---|---|
| Feeding system and corresponding activity coefficients | ||||
| Feeding system (1 Stall - 2 Pasture - 3 Grazing large area) | 1.00 | 1.00 | const | feeding_system |
| Acticity coefficient Ca - Stall feeding | 0.00 | 0.00 | const | Ca_stall |
| Acticity coefficient Ca -Pasture feeding | 0.14 | 0.20 | posnorm | Ca_pasture |
| Acticity coefficient Ca - Grazing large area | 0.31 | 0.41 | posnorm | Ca_large_grazing |
| Herd composition | ||||
| Number of ‘intact’ mature bulls (>36 months) (animal type 1) | 5.00 | 5.00 | const | N_animal_type_1 |
| Number of mature castrates (>36 months) (animal type 2) | 2.00 | 2.00 | const | N_animal_type_2 |
| Number of males 12-36 months (animal type 3) | 4.00 | 4.00 | const | N_animal_type_3 |
| Number of lactating cows (animal type 4) | 5.00 | 5.00 | const | N_animal_type_4 |
| Number of lactating cows (animal type 5) | 3.00 | 3.00 | const | N_animal_type_5 |
| Number of lactating cows (animal type 6) | 2.00 | 2.00 | const | N_animal_type_6 |
| Number of weaned young females <=12 (animal type 7) | 4.00 | 4.00 | const | N_animal_type_7 |
| Number of weaned young males <=12 (animal type 8) | 3.00 | 3.00 | const | N_animal_type_8 |
| Number of females 12-36 months (animal type 9) | 5.00 | 5.00 | const | N_animal_type_9 |
| Number of pre-weaned female calves (animal type 10) | 6.00 | 6.00 | const | N_animal_type_10 |
| Number of pre-weaned male calves (animal type 10) | 4.00 | 4.00 | const | N_animal_type_11 |
| Number of first-time pregnant cows (animal type 12) | 3.00 | 3.00 | const | N_animal_type_12 |
| Sex-specific and (maintenance) coefficients | ||||
| Sex-specific coefficient - intact males | 1.02 | 1.38 | posnorm | C_intact_males |
| Sex-specific coefficient - castrates | 0.85 | 1.15 | posnorm | C_castrates |
| Sex-specific coefficient - females | 0.68 | 0.92 | posnorm | C_females |
| Maintenance coefficient (Cfi) - Cattle/Buffalo (non-lactating cows, steers and juveniles) | 0.27 | 0.37 | posnorm | maintenance_coefficient_Cfi_other |
| Maintenance coefficient (Cfi) - Cattle/Buffalo (lactating cows) | 0.33 | 0.44 | posnorm | maintenance_coefficient_Cfi_lact_cow |
| Maintenance coefficient (Cfi) - Cattle/Buffalo (bulls) | 0.31 | 0.43 | posnorm | maintenance_coefficient_Cfi_bulls |
| Mature body weight females | 279.20 | 418.80 | posnorm | mature_weight_MW_female |
| Mature body weight males | 386.40 | 579.60 | posnorm | mature_weight_MW_male |
| Weight (gain) data by animal type | ||||
| Mean live weight of animal type 1 | 99.50 | 537.68 | posnorm | live_weight_LW_1 |
| Mean live weight of animal type 2 | 84.40 | 468.84 | posnorm | live_weight_LW_2 |
| Mean live weight of animal type 3 | 81.50 | 374.07 | posnorm | live_weight_LW_3 |
| Mean live weight of animal type 4 | 69.20 | 480.83 | posnorm | live_weight_LW_4 |
| Mean live weight of animal type 5 | 59.20 | 454.38 | posnorm | live_weight_LW_5 |
| Mean live weight of animal type 6 | 62.40 | 488.65 | posnorm | live_weight_LW_6 |
| Mean live weight of animal type 7 | 61.90 | 233.83 | posnorm | live_weight_LW_7 |
| Mean live weight of animal type 8 | 45.20 | 191.35 | posnorm | live_weight_LW_8 |
| Mean live weight of animal type 9 | 90.80 | 402.37 | posnorm | live_weight_LW_9 |
| Mean live weight of animal type 10 | 23.60 | 93.12 | posnorm | live_weight_LW_10 |
| Mean live weight of animal type 11 | 18.90 | 96.89 | posnorm | live_weight_LW_11 |
| Mean live weight of animal type 12 | 63.90 | 425.82 | posnorm | live_weight_LW_12 |
| Average daily weight gain (kg/day) of animal type 1 | 0.00 | 0.00 | const | weight_gain_WG_1 |
| Average daily weight gain (kg/day) of animal type 2 | 0.00 | 0.00 | const | weight_gain_WG_2 |
| Average daily weight gain (kg/day) of animal type 3 | 0.12 | 0.18 | posnorm | weight_gain_WG_3 |
| Average daily weight gain (kg/day) of animal type 4 | 0.00 | 0.00 | const | weight_gain_WG_4 |
| Average daily weight gain (kg/day) of animal type 5 | 0.00 | 0.00 | const | weight_gain_WG_5 |
| Average daily weight gain (kg/day) of animal type 6 | 0.00 | 0.00 | const | weight_gain_WG_6 |
| Average daily weight gain (kg/day) of animal type 7 | 0.34 | 0.50 | posnorm | weight_gain_WG_7 |
| Average daily weight gain (kg/day) of animal type 8 | 0.28 | 0.40 | posnorm | weight_gain_WG_8 |
| Average daily weight gain (kg/day) of animal type 9 | 0.21 | 0.31 | posnorm | weight_gain_WG_9 |
| Average daily weight gain (kg/day) of animal type 10 | 0.34 | 0.50 | posnorm | weight_gain_WG_10 |
| Average daily weight gain (kg/day) of animal type 11 | 0.28 | 0.40 | posnorm | weight_gain_WG_11 |
| Average daily weight gain (kg/day) of animal type 12 | 0.21 | 0.31 | posnorm | weight_gain_WG_12 |
| Milk yield; fat content; pregnancy coefficients | ||||
| Average milk yield (kg/day) of animal type 4 | 5.11 | 8.35 | posnorm | milk_yield_4 |
| Average milk yield (kg/day) of animal type 5 | 5.11 | 8.35 | posnorm | milk_yield_5 |
| Average milk yield (kg/day) of animal type 6 | 5.11 | 8.35 | posnorm | milk_yield_6 |
| Milk fat (%) | 3.60 | 4.40 | posnorm | milk_fat |
| Maximum pregnancy coefficient Cp (scaled to year) | 0.10 | 0.10 | const | Cp_0 |
| Pregnancy coefficient for animal type 4 (share of full coefficient) | 0.00 | 0.89 | posnorm | Cp_rate_4 |
| Pregnancy coefficient for animal type 5 (share of full coefficient) | 0.11 | 1.00 | posnorm | Cp_rate_5 |
| Pregnancy coefficient for animal type 6 (share of full coefficient) | 0.01 | 0.99 | posnorm | Cp_rate_6 |
| Emissions from externally produced feed | ||||
| Feed production | 834.56 | 834.56 | const | feed_prod_CO2 |
| Feed transport | 5.20 | 5.20 | const | feed_trans_CO2 |
| Conversion parameters | ||||
| Digestible energy expressed as a percentage of gross energyy (De) | 58.67 | 59.04 | posnorm | DE_perc |
| Percent crude protein in diet | 5.00 | 15.00 | posnorm | perc_crude_protein_diet |
| Methane conversion rate (0-1) | 0.05 | 0.08 | posnorm | Y_m |
| emission factor for N2O emissions from nitrogen leaching and runoff, kg N2O-N/kg N leached and runoff | 0.01 | 0.01 | const | EF5 |
| emission factor for N2O emissions from atmospheric deposition of nitrogen on soils and water surfaces | 0.01 | 0.01 | const | EF4 |
| urinary energy expressed as fraction of Gross Energy intake | 0.04 | 0.04 | const | urinary_energy |
| the ash content of manure calculated as a fraction of the dry matter feed intake | 0.08 | 0.08 | const | ash_content |
| Constants for methane production and conversion | ||||
| Maximum methane producing capacity of the manure (animal type 1) | 0.11 | 0.15 | posnorm | Bo_animaltype_1 |
| Maximum methane producing capacity of the manure (animal type 2) | 0.11 | 0.15 | posnorm | Bo_animaltype_2 |
| Maximum methane producing capacity of the manure (animal type 3) | 0.11 | 0.15 | posnorm | Bo_animaltype_3 |
| Maximum methane producing capacity of the manure (animal type 4) | 0.11 | 0.15 | posnorm | Bo_animaltype_4 |
| Maximum methane producing capacity of the manure (animal type 5) | 0.11 | 0.15 | posnorm | Bo_animaltype_5 |
| Maximum methane producing capacity of the manure (animal type 6) | 0.11 | 0.15 | posnorm | Bo_animaltype_6 |
| Maximum methane producing capacity of the manure (animal type 7) | 0.11 | 0.15 | posnorm | Bo_animaltype_7 |
| Maximum methane producing capacity of the manure (animal type 8) | 0.11 | 0.15 | posnorm | Bo_animaltype_8 |
| Maximum methane producing capacity of the manure (animal type 9) | 0.11 | 0.15 | posnorm | Bo_animaltype_9 |
| Maximum methane producing capacity of the manure (animal type 10) | 0.11 | 0.15 | posnorm | Bo_animaltype_10 |
| Maximum methane producing capacity of the manure (animal type 11) | 0.11 | 0.15 | posnorm | Bo_animaltype_11 |
| Maximum methane producing capacity of the manure (animal type 12) | 0.11 | 0.15 | posnorm | Bo_animaltype_12 |
| Methane conversion factor - pasture | 0.00 | 0.01 | posnorm | MCFpasture |
| Methane conversion factor - spreading | 0.00 | 0.01 | posnorm | MCFspread |
| Methane conversion factor - drylot | 0.01 | 0.03 | posnorm | MCFdrylot |
| Methane conversion factor - solid storage | 0.03 | 0.07 | posnorm | MCFsolidstore |
| Methane conversion factor - composted | 0.01 | 0.04 | posnorm | MCFcomposted |
| Methane conversion factor - liquid | 0.40 | 1.20 | posnorm | MCFliquid |
| Methane conversion factor - biogas | 0.05 | 0.16 | posnorm | MCFbiogas |
| Methane conversion factor - burning | 0.05 | 0.15 | posnorm | MCFburn |
| Methane conversion factor - sold | 0.00 | 0.00 | const | MCFsold |
| Decision tree for manure management (no longer needed) | ||||
| question 1: how much of the manure is centrally collected? | 0.10 | 0.90 | tnorm_0_1 | share_centrally |
| question 2: of manure that’s not centrally collected, how much is burned? | 0.10 | 0.90 | tnorm_0_1 | share_central_burn |
| question 3: of centrally collected manure, how much is spread daily on cropland or pasture? | 0.10 | 0.90 | tnorm_0_1 | share_spread |
| question 4: how much of the manure that’s not spread is sold? | 0.10 | 0.90 | tnorm_0_1 | share_sold |
| question 5: how much of the stored manure that’s not sold is stored in liquid form? | 0.10 | 0.90 | tnorm_0_1 | share_liquid |
| question 6: how much of the liquid manure is converted to biogas? | 0.10 | 0.90 | tnorm_0_1 | share_biogas |
| question 7: how much of the solid manure is composted? | 0.10 | 0.90 | tnorm_0_1 | share_compost |
| question 8: how much of the solid manure is stored in dry lots? | 0.10 | 0.90 | tnorm_0_1 | share_drylot |
| Manure management system | ||||
| Share of manure management (pasture) | 0.30 | 0.30 | const | mms_pasture |
| Share of manure management (spread) | 0.10 | 0.10 | const | mms_spread |
| Share of manure management (dry lot) | 0.10 | 0.10 | const | mms_drylot |
| Share of manure management (solid storage) | 0.10 | 0.10 | const | mms_solidstore |
| Share of manure management (composted) | 0.10 | 0.10 | const | mms_composted |
| Share of manure management (liquid) | 0.10 | 0.10 | const | mms_liquid |
| Share of manure management (biogas) | 0.10 | 0.10 | const | mms_biogas |
| Share of manure management (burned) | 0.10 | 0.10 | const | mms_burn |
| Share of manure management (burned) | 0.00 | 0.00 | const | mms_sold |
| Manure emission coefficients by management | ||||
| Manure nitrogen direct emission coefficients (Table 10.21 from the IPCC guidelines) | ||||
| Manure nitrogen direct emission coefficient (pasture) | 0.00 | 0.04 | posnorm | mndec_pasture |
| Manure nitrogen direct emission coefficient (spread) | 0.00 | 0.00 | const | mndec_spread |
| Manure nitrogen direct emission coefficient (dry lot) | 0.00 | 0.04 | posnorm | mndec_drylot |
| Manure nitrogen direct emission coefficient (solid storage) | 0.00 | 0.02 | posnorm | mndec_solidstore |
| Manure nitrogen direct emission coefficient (composted) | 0.00 | 0.02 | posnorm | mndec_composted |
| Manure nitrogen direct emission coefficient (liquid) | 0.00 | 0.02 | posnorm | mndec_liquid |
| Manure nitrogen direct emission coefficient (biogas) | 0.00 | 0.00 | const | mndec_biogas |
| Manure nitrogen direct emission coefficient (burned) | 0.00 | 0.00 | const | mndec_burn |
| Manure nitrogen direct emission coefficient (sold) | 0.00 | 0.00 | const | mndec_sold |
| Manure nitrogen volatilization emission coefficients (Table 10.22 from the IPCC guidelines) | ||||
| Manure nitrogen volatilization emission coefficient (pasture) | 0.02 | 0.02 | const | mnvc_pasture |
| Manure nitrogen volatilization emission coefficient (spread) | 0.05 | 0.60 | posnorm | mnvc_spread |
| Manure nitrogen volatilization emission coefficient (dry lot) | 0.20 | 0.50 | posnorm | mnvc_drylot |
| Manure nitrogen volatilization emission coefficient (solid storage) | 0.10 | 0.65 | posnorm | mnvc_solidstore |
| Manure nitrogen volatilization emission coefficient (composted) | 0.14 | 0.70 | posnorm | mnvc_composted |
| Manure nitrogen volatilization emission coefficient (liquid) | 0.09 | 0.36 | posnorm | mnvc_liquid |
| Manure nitrogen volatilization emission coefficient (biogas) | 0.00 | 0.00 | const | mnvc_biogas |
| Manure nitrogen volatilization emission coefficient (burned) | 0.00 | 0.00 | const | mnvc_burn |
| Manure nitrogen volatilization emission coefficient (sold) | 0.00 | 0.00 | const | mnvc_sold |
| Fraction of managed manure nitrogen losses for livestock category T due to runoff and leaching during solid and liquid storage of manure (typical range 0.01-0.20) | ||||
| Fraction of managed manure nitrogen losses … (pasture) | 0.01 | 0.73 | posnorm | frac_leach_pasture |
| Fraction of managed manure nitrogen losses … (spread) | 0.01 | 0.73 | posnorm | frac_leach_spread |
| Fraction of managed manure nitrogen losses … (dry lot) | 0.01 | 0.73 | posnorm | frac_leach_drylot |
| Fraction of managed manure nitrogen losses … (solid storage) | 0.01 | 0.73 | posnorm | frac_leach_solidstore |
| Fraction of managed manure nitrogen losses … (composted) | 0.01 | 0.73 | posnorm | frac_leach_composted |
| Fraction of managed manure nitrogen losses … (liquid) | 0.01 | 0.73 | posnorm | frac_leach_liquid |
| Fraction of managed manure nitrogen losses … (biogas) | 0.00 | 0.00 | const | frac_leach_biogas |
| Fraction of managed manure nitrogen losses … (burned) | 0.00 | 0.00 | const | frac_leach_burn |
| Fraction of managed manure nitrogen losses … (sold) | 0.00 | 0.00 | const | frac_leach_sold |
Monte Carlo simulation
With the model equation and the input table defined above, we can run a Monte Carlo simulation that returns a probability distribution for the greenhouse gas emissions resulting from the particular herd that is specified in the input table. We use the functions provided by the decisionSupport package (Luedeling et al. 2022) in the R programming language (R Core Team 2022).
GHG_emissions<-mcSimulation(estimate=estimate_read_csv(input_table),
model_function=ghg_emissions,
numberOfModelRuns=1000,
functionSyntax="plainNames")Running a Monte Carlo simulation with scenarios
The mcSimulation function takes random draws for each of the input variables from the distributions specified in the input table. The distributions that are used are the same for every run of the Monte Carlo simulation. This may be desirable in homogeneous populations where all combinations of random values are plausible, but it often misses the mark where populations are composed of particular member types with specific characteristics. To account for populations that are heterogeneous, which is the case for livestock-keeping households in Kenya, we use the scenario_mc function, which allows specifying particular subgroups. Defining such scenarios can be desirable when it’s necessary to simulate outcomes for particular farm types, climate scenarios etc. To do this, we would normally have to define separate input tables and run separate simulations. This is quite inconvenient. Here, we use a function that can do this job automatically.
The scenario_mc function in the decisionSupport package (Luedeling et al. 2022) consists of a wrapper around the mcSimulation function. The function is able to modify the input table according to a scenario file we provide. This file contains information on which variables should be modified for each scenario and how many model runs should be included in the simulation. The function is essentially the same as the mcSimulation function, but with an additional scenarios option which imports the contents of a .csv file.
The .csv file for the scenarios option specifies variables that should be modified for each of the scenarios, as well as the number of model runs for it. Each column beyond the first two columns in the file contain information for a particular scenario. This is restricted to the parameters that should be changed for this scenario.
Here’s an example scenarios file:
knitr::kable(read.table("data/scenarios.csv", sep = ",", fileEncoding = "UTF-8-BOM",header=TRUE))| Variable | param | Scen1 | Scen2 |
|---|---|---|---|
| Runs | x | 100 | 300 |
| N_animal_type_1 | both | 100 | 4000 |
| N_animal_type_2 | lower | 100 | 10 |
| N_animal_type_2 | upper | 1000 | 20 |
| N_animal_type_2 | distribution | posnorm | norm |
The first row (“Runs” in the first column) indicates the number of runs. Each subsequent row indicates which value from the input table should be modified, so the string in the first column has to correspond to a variable name in the input table. The second column indicates which values should be changed, with options being “lower,” “upper,” “distribution” and “both.” If “both” is selected, both the lower and upper bounds are modified (this should only be done for variables witha constant distribution, since the rplacement will be for both the upper and lower bound).
We define subgroups based on the populations of households described by a household survey in Kenya (Wilkes et al. 2020, 2019). Each of the farms encountered in this survey were converted into one scenario, for which Monte Carlo simulation is then used to produce multiple replicates.
Household-based scenarios
The information for the household is extracted from a spreadsheet that contains all the survey responses (Wilkes et al. 2020, 2019). Extracting all relevant information from these spreadsheets is a slightly complex task, which is accomplished by the following code.
# import the data from the survey and ensure that the .csv encoding can be read by R
Kenyaherds<-read.table("data/Kenya_baseline_individual_cow_wNotes.csv", sep = ",",
fileEncoding = "UTF-8-BOM", header=TRUE)
# count the number of animals for each type
animals<-table(Kenyaherds[,c("QQNO","animal_type")])
variables_to_update_const<-c(paste0("N_animal_type_",1:12),
"feeding_system")
scenarios<-data.frame(Variable=variables_to_update_const,param="both")
variables_to_update_posnorm<-c(paste0("milk_yield_",4:6),
paste0("Cp_rate_",4:6),
"feed_prod_CO2",
"feed_trans_CO2",
paste0("live_weight_LW_",1:12),
paste0("weight_gain_WG_",1:12),
"DE_perc",
"perc_crude_protein_diet")
scenarios<-rbind(scenarios,data.frame(Variable=rep(variables_to_update_posnorm,each=3),param=c("lower","upper","distribution")))
for(i in 1:nrow(animals))
{coln<-paste0("Farm_",row.names(animals)[i])
scenarios[,coln]<-NA
scenarios[which(scenarios$Variable %in% paste0("N_animal_type_",1:12)),coln]<-animals[i,]
}
# for some animal-scale variables (e.g. milk yield), we need averages per farm
### error estimates for variables
milk_yield_error<-0.275 # Migose et al. (2020, https://www.sciencedirect.com/science/article/abs/pii/S1871141319307371 )
feed_trans_error<-0.25 # Vellinga et al (2013, https://edepot.wur.nl/254098)
feed_prod_error<-0.25 # Vellinga et al (2013, https://edepot.wur.nl/254098)
live_weight_error<-0.145 # Goopy et al. (2018, https://www.publish.csiro.au/an/an16577 )
weight_gain_error<-0.145 # Goopy et al. (2018, https://www.publish.csiro.au/an/an16577 )
digestible_energy_error<-0.1
crude_protein_error<-0.1
# farm-scale variables
Farms<-unique(Kenyaherds$QQNO)
for (f in Farms)
{coln<-paste0("Farm_",f)
# update feeding system
scenarios[which(scenarios$Variable == "feeding_system"),coln]<-
median(Kenyaherds[which(Kenyaherds$QQNO==f),"system"])
# update digestible energy (mean for all animals for now) (with error)
scenarios[which(scenarios$Variable == "DE_perc"),coln][1:2]<-
mean(Kenyaherds[which(Kenyaherds$QQNO==f),"dig_energy"])*
c(1-digestible_energy_error,1+digestible_energy_error)
# update crude protein (mean for all animals for now) (with error)
scenarios[which(scenarios$Variable == "perc_crude_protein_diet"),coln][1:2]<-
mean(Kenyaherds[which(Kenyaherds$QQNO==f),"CP."])*
c(1-crude_protein_error,1+crude_protein_error)
# update feed production CO2 (with error)
scenarios[which(scenarios$Variable == "feed_prod_CO2"),coln][1:2]<-
mean(Kenyaherds[which(Kenyaherds$QQNO==f),"feed_kgCO2"])*
c(1-feed_prod_error,1+feed_prod_error)
# update feed transport CO2 (with error)
scenarios[which(scenarios$Variable == "feed_trans_CO2"),coln][1:2]<-
mean(Kenyaherds[which(Kenyaherds$QQNO==f),"feed.transport_kgCO2"])*
c(1-feed_trans_error,1+feed_trans_error)
# update milk yield for each cow type (animal type 4-6) (with error)
for(typ in 4:6)
if(sum(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ)>0)
{scenarios[which(scenarios$Variable==paste0("milk_yield_",typ)),coln][1:2]<-
mean(Kenyaherds[which(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ),
"milk_yield"])*c(1-milk_yield_error,1+milk_yield_error)
scenarios[which(scenarios$Variable==paste0("CP_rate_",typ)),coln][1:2]<-
mean(Kenyaherds[which(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ),
"Cp"]/0.1) *c(1,1)
}
# update live weight and weight gain for each animal type (with error)
for(typ in 1:12)
if(sum(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ)>0)
{scenarios[which(scenarios$Variable==paste0("live_weight_LW_",typ)),coln][1:2]<-
mean(Kenyaherds[which(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ),
"live_weight_LW"])*c(1-live_weight_error,1+live_weight_error)
scenarios[which(scenarios$Variable==paste0("weight_gain_WG_",typ)),coln][1:2]<-
mean(Kenyaherds[which(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ),
"weight_gain_WG"])*c(1-weight_gain_error,1+weight_gain_error)
}
for(varia in variables_to_update_posnorm)
if(sum(is.na(scenarios[which(scenarios$Variable==varia),coln][1:2]))==0)
if(!equals(scenarios[which(scenarios$Variable==varia),coln][1],
scenarios[which(scenarios$Variable==varia),coln][2]))
scenarios[which(scenarios$Variable==varia&
scenarios$param=="distribution"),coln]<-"posnorm" else
scenarios[which(scenarios$Variable==varia&
scenarios$param=="distribution"),coln]<-"const"
}
write.csv(scenarios,"data/farm_scenarios.csv",row.names = FALSE)Here’s the resulting file (only the first 3 farms are shown, and only the first 20 rows):
knitr::kable(read.table("data/farm_scenarios.csv", sep = ",", fileEncoding = "UTF-8-BOM",header=TRUE)[1:20,1:5])| Variable | param | Farm_1 | Farm_2 | Farm_3 |
|---|---|---|---|---|
| N_animal_type_1 | both | 0 | 1 | 1 |
| N_animal_type_2 | both | 0 | 0 | 0 |
| N_animal_type_3 | both | 0 | 0 | 0 |
| N_animal_type_4 | both | 0 | 1 | 0 |
| N_animal_type_5 | both | 0 | 0 | 2 |
| N_animal_type_6 | both | 2 | 0 | 0 |
| N_animal_type_7 | both | 1 | 0 | 0 |
| N_animal_type_8 | both | 0 | 0 | 0 |
| N_animal_type_9 | both | 0 | 1 | 0 |
| N_animal_type_10 | both | 0 | 0 | 0 |
| N_animal_type_11 | both | 0 | 0 | 0 |
| N_animal_type_12 | both | 0 | 0 | 0 |
| feeding_system | both | 2 | 3 | 3 |
| milk_yield_4 | lower | 2.826069862825 | ||
| milk_yield_4 | upper | 4.969984931175 | ||
| milk_yield_4 | distribution | posnorm | ||
| milk_yield_5 | lower | 3.8373469091625 | ||
| milk_yield_5 | upper | 6.7484376678375 | ||
| milk_yield_5 | distribution | posnorm | ||
| milk_yield_6 | lower | 5.471348045275 |
The first few rows in this table adjust the herd composition parameters, so that the scenario herd corresponds to the herd of the specific household. Other parameters including the milk yield per cow are also adjusted.
Now we can run the Monte Carlo simulation with these scenarios.
GHG_simulation_scenarios<-scenario_mc(base_estimate = estimate_read_csv("data/livestock_ghg_input_table.csv"),
scenarios = read.csv("data/farm_scenarios.csv",
fileEncoding = "UTF-8-BOM"),
model_function = ghg_emissions,
numberOfModelRuns = 10,
functionSyntax = "plainNames")Given these results we can now plot distributions of the various outputs from the model, using the plot_distributions function.
# enteric_CH4=enteric_CH4,
enteric_CH4_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
vars = "enteric_CH4",
method = "boxplot_density",
old_names = "enteric_CH4",
new_names = "CH4 Emissions from Enteric Fermentation per Farm") +
theme(axis.title.x=element_blank(),
axis.title.y=element_blank())
# milk_yield=sum(total_milk),
milk_yield_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
vars = "milk_yield",
method = "boxplot_density",
old_names = "milk_yield",
new_names = "Total Milk Yield") +
theme(axis.title.x=element_blank(),
axis.title.y=element_blank())
# enteric_CH4_CO2eq=sum(enteric_CH4_CO2eq),
enteric_CH4_CO2eq_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
vars = "enteric_CH4_CO2eq",
method = "boxplot_density",
old_names = "enteric_CH4_CO2eq",
new_names = "Carbon Dioxide Equivalent Enteric Methane Emissions from Livestock") +
theme(axis.title.x=element_blank())
# mm_CH4_CO2eq=sum(mm_CH4_CO2eq),
mm_CH4_CO2eq_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
vars = "mm_CH4_CO2eq",
method = "boxplot_density",
old_names = "mm_CH4_CO2eq",
new_names = "Carbon Dioxide Equivalent Methane Emissions from Livestock Manure") +
theme(axis.title.x=element_blank(),
axis.title.y=element_blank())
# mm_N2O_CO2eq=sum(mm_N2O_CO2eq),
mm_N2O_CO2eq_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
vars = "mm_N2O_CO2eq",
method = "boxplot_density",
old_names = "mm_N2O_CO2eq",
new_names = "Carbon Dioxide Equivalent Nitrous Oxide Emissions from Livestock Manure") +
theme(axis.title.x=element_blank(),
axis.title.y=element_blank())
# feed_CO2=feed_CO2,
feed_CO2_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
vars = "feed_CO2",
method = "boxplot_density",
old_names = "feed_CO2",
new_names = "Carbon Dioxide Emissions from Livestock Feed Production and Transport") +
theme(axis.title.x=element_blank(),
axis.title.y=element_blank())
# on_farm=sum(on_farm)
on_farm_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
vars = "on_farm",
method = "boxplot_density",
old_names = "on_farm",
new_names = "Total Carbon Dioxide Equivalent Emissions from Livestock") +
theme(axis.title.y=element_blank())
library(patchwork)##
## Attaching package: 'patchwork'## The following object is masked from 'package:raster':
##
## area# remove the axis marks from the two lefthand-most plots and the x-axis title from the left and right plots
# plot with patchwork
layout <- "
AAABBB
CCCDDD
EEEFFF
GGGGGG
"
enteric_CH4_dist_plot + milk_yield_dist_plot +
enteric_CH4_CO2eq_dist_plot + mm_CH4_CO2eq_dist_plot +
mm_N2O_CO2eq_dist_plot + feed_CO2_dist_plot +
on_farm_dist_plot +
plot_layout(guides = "collect", design = layout) 
ggsave("Fig_1_GHG_simulation_outputs.png")## Saving 10 x 10 in image# THIS STILL NEEDS IMPROVEMENT - e.g. units are missing, and plot titles don’t fit (use concise symbols with subscripts etc.)
Possible errors resulting from input uncertainty
To illustrate the possible influence of input uncertainty on estimates of total farm emissions, we show the example of farm #5 from the household survey. This farm has 6 head of cattle, including 3 lactating cows. We simulate 1000 versions of this farm.
scenarios<-read.csv("data/farm_scenarios.csv", fileEncoding = "UTF-8-BOM")
farm_5<-scenarios[,c(1,2,7)]
farm5_emissions<-scenario_mc(base_estimate = estimate_read_csv("data/livestock_ghg_input_table.csv"),
scenarios = farm_5,
model_function = ghg_emissions,
numberOfModelRuns = 1000,
functionSyntax = "plainNames")
plot_distributions(mcSimulation_object = farm5_emissions,
vars = "on_farm",
method = "boxplot_density",
old_names = "on_farm",
new_names = "Total Carbon Dioxide Equivalent Emissions from Livestock") +
theme(axis.title.y=element_blank())
ggsave("Fig_2_uncertainty_illustration.png")## Saving 7 x 5 in image# NEED TO ADD THE USUAL PRECISE VALUE HERE FOR COMPARISON (also need to add units and possibly fix up other things)
This
farm5_emissions$y[,"CO2eq_per_milk"]<-farm5_emissions$y$pc_on_farm/farm5_emissions$y$pc_milk_yield
ggplot(data=farm5_emissions$y, aes(CO2eq_per_milk)) +
geom_histogram(aes(y=..density..), colour="black", fill="white")+
geom_density(alpha=.2, fill="#FF6666") +
ylab("Relative probability") +
xlab(expression(Emission~density~of~milk~production~(kg~CO[2]-eq~kg^-1)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggsave("Fig_3_emission_density_farm5.png")## Saving 7 x 5 in image
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.# STILL NEED TO add the usual precise value estimate
Population-wide emission density of milk production
Since the emission density of each household thus involves considerable uncertainty, such uncertainty should also be expected at the population level.
We can illustrate this uncertainty using a scatter plot. The following figure shows all combinations of total milk yield and total greenhouse gas emissions across all households.
As with the GHG Model, we apply the uncertainty functions on the milk_yield and on_farm variables. These come from the model results (4140) of the GHG_simulation_scenarios model results to the in_var (x) and an outcome variable out_var (y). In the model we defined on_farm as the sum of Carbon Dioxide Equivalent Enteric Methane Emissions from Livestock (enteric_CH4_CO2eq), Carbon Dioxide Equivalent Methane Emissions from Livestock Manure (mm_CH4_CO2eq), Carbon Dioxide Equivalent Nitrous Oxide Emissions from Livestock Manure (mm_N2O_CO2eq) and Carbon Dioxide Emissions from Livestock Feed Production and Transport (feed_CO2).
uncertainty::varscatter(in_var = GHG_simulation_scenarios$y$pc_milk_yield,
out_var = GHG_simulation_scenarios$y$pc_on_farm,
xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1})))
## [1] "Scatter plot of influencing variable (x) and outcome variable (y) with associated and estimated densities."# NICE PLOT, but would be nicer still if this was also a ggplot figure. Kinda inconsistent if it isn’t. Also can’t save it with ggsave, so will need another way.
Emissions intensity of milk production
An indicator that is often of interest is the ratio between these two variables: total emissions divided by the milk yield. We’ll calculate this next. We then have to remove all households that don’t produce any milk, because this ratio would be infinity for them. In fact, we’ll remove all households that reported very low milk yields of <500 kg per year per capita, since these are unlikely to represent regularly operating dairy operations. For instance, they may be mostly focused on beef cattle production, or they may only have a small number of dairy cows that happened to be unproductive at the time of the survey.
GHG_simulation_scenarios$y[,"CO2eq_per_milk"]<-GHG_simulation_scenarios$y$pc_on_farm/GHG_simulation_scenarios$y$pc_milk_yield
normal_milk_yield<-which(GHG_simulation_scenarios$y$pc_milk_yield>500)
dairy_households<-GHG_simulation_scenarios
dairy_households$x<-dairy_households$x[normal_milk_yield,]
dairy_households$y<-dairy_households$y[normal_milk_yield,]Now we can plot the emissions intensity vs. the milk yield.
uncertainty::varscatter(in_var = dairy_households$y$pc_milk_yield,
out_var = dairy_households$y$CO2eq_per_milk,
xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
ylab = expression(Total~emissions~intensity~of~milk~production~(kg~CO[2]-eq~kg^{-1})))
## [1] "Scatter plot of influencing variable (x) and outcome variable (y) with associated and estimated densities."This can also be illustrated by a density surface.
uncertainty::varkernel(in_var = dairy_households$y$pc_milk_yield,
out_var = dairy_households$y$CO2eq_per_milk,
xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
ylab = expression(Total~emissions~intensity~of~milk~production~(kg~CO[2]-eq~kg^{-1})))
## [1] "Density surface plot of an expected influential variable (x) and an outcome variable of interest (y)."The milk yield per capita is sometimes proposed as an indicator of greenhouse gas emissions. We can now examine the merits of this indicator. Essentially, its usefulness hinges on its relationship with greenhouse gas emission intensity. We can investigate this by taking slices through the density surface.
Probability of GHG emissions for a given milk yield
The following plot shows the emissions intensity that corresponds to a milk yield of 2000 kg per capita per year, based on the population of households covered in the survey.
slice<-uncertainty::varslice_resample(in_var = dairy_households$y$pc_milk_yield,
out_var = dairy_households$y$CO2eq_per_milk,
expectedin_var=2000,
n = 1000,
n_samples = 1000,
out_var_sampling = 1000)
graphics::plot(x = slice$slice$Output_values,
y = slice$slice$Relative_probability,
type = "l",
ylab = "Relative probability",
xlab = "Emissions intensity",
col = "seagreen",
lwd = 2)
ggplot(data=data.frame(CO2eq_per_milk=slice$resampled), aes(CO2eq_per_milk)) +
geom_histogram(aes(y=..density..), colour="black", fill="white")+
geom_density(alpha=.2, fill="#FF6666") +
ylab("Relative probability") +
xlab(expression(Emission~density~of~milk~production~at~2000~kg~milk~per~capita~(XX~CO[2]-eq~kg^-1)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggsave("Fig_X_emission_density_all_farms.png")## Saving 7 x 5 in image
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.We see considerable uncertainty about the emissions at this particular level of production, and the previous figure shows that emissions of all other levels is similarly variable. This raises considerable doubts about the usefulness of milk yield as an indicator of emissions intensity - at least when this indicator is taken at face value, without considering this uncertainty. Yet, some information is contained even in such an imperfect indicator, since, as also shown in the previous figure, there is clearly a correlation between milk yield and emissions intensity. The following proposes a mechanism for quantifying the level of certainty, with which a change in milk yield can reasonably be associated with a change in emissions intensity.
Probability of an increase in milk yield being associated with an increase in GHG emission intensity
To derive the level of certainty, with which a change in milk yield is associated with an increase or decrease in emissions intensity, we have to investigate the probability distributions (as shown in Fig. X) of emissions intensity at two levels of milk production. This can be achieved using the varslice_resample function for both production levels.
Let’s look at this for production of 2000 and 3000 kg milk per capita per year. We’ll take 1000 random draws of the emissions intensities that associated with both levels and compare values among these samples.
slice_2000<-uncertainty::varslice_resample(in_var = dairy_households$y$pc_milk_yield,
out_var = dairy_households$y$CO2eq_per_milk,
expectedin_var=2000,
n = 1000,
n_samples = 1000,
out_var_sampling = 1000)
slice_3000<-uncertainty::varslice_resample(in_var = dairy_households$y$pc_milk_yield,
out_var = dairy_households$y$CO2eq_per_milk,
expectedin_var=3000,
n = 1000,
n_samples = 1000,
out_var_sampling = 1000)
sum(slice_3000$resampled < slice_2000$resampled) / length(slice_2000$resampled)## [1] 0The comparison in the last line of code above reveals that in 0 % of cases, emissions intensity at the higher milk production level is lower than at the lower level.
We can now generalize this result by doing similar calculations for multiple production levels and illustrating them for the whole population.
in_var = dairy_households$y$pc_milk_yield
out_var = dairy_households$y$CO2eq_per_milk
milk_yields_to_sample<-seq(100,10000,100)
milk_yields_to_sample<-milk_yields_to_sample[which(milk_yields_to_sample<max(in_var))]
milk_yields_to_sample<-milk_yields_to_sample[which(milk_yields_to_sample>min(in_var))]
slices<-list()
# slice through the distribution at every 100 L milk yield per cow
for(i in 1:length(milk_yields_to_sample))
{slices[[i]]<-
varslice_resample(in_var,
out_var,
expectedin_var = milk_yields_to_sample[i],
n_samples = 100000,
out_var_sampling = 1000)$resample
}
decreased_emissions_probability <-
data.frame(Base_yield = NA,new_yield=NA,probability=NA)
for (i in 1:(length(milk_yields_to_sample) - 1))
{
for (j in (i + 1):length(milk_yields_to_sample))
{decreased_emissions_probability[nrow(decreased_emissions_probability)+1,"Base_yield"]<-milk_yields_to_sample[i]
decreased_emissions_probability$new_yield[nrow(decreased_emissions_probability)]<-milk_yields_to_sample[j]
decreased_emissions_probability$probability[nrow(decreased_emissions_probability)] <-
sum(slices[[i]] > slices[[j]]) / length(slices[[i]])
}
}
decreased_emissions_probability<-decreased_emissions_probability[which(!is.na(decreased_emissions_probability$Base_yield)),]
decreased_emissions_probability[,"yield_increase"]<-decreased_emissions_probability$new_yield-decreased_emissions_probability$Base_yield
ggplot(data = decreased_emissions_probability, aes(Base_yield, yield_increase, z =
probability)) + geom_contour_filled() +
ggtitle("Probability that emission intensity reductions occurred") +
xlab("Baseline milk yield per head and year (kg)") +
ylab("Milk yield increase per head and year (kg)") +
stat_contour(breaks = c(0.9),col="black")## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf
ggsave("Fig_X1_probability_of_intensity_reduction.png")## Saving 7 x 5 in image## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -InfThe 90% confidence contour is shown in the plot. The results indicate that whenever milk yields increase, the probability that emissions intensity decreased is usually greater than 50% (slightly lower values are due to randomness in the resampling procedure), but that large increases are needed to be sure of this.
Testing interventions - strategy
Farms may try to reduce their emissions intensity by making changes to their herd structure, feeding strategy or management practices. Programs aiming to effect emissions reductions in dairy production may want to reward such changes, and they may be looking to the milk production per capita as an indicator for the effectiveness of such measures. We now want explore whether this is a reasonable strategy.
Strategy
Develop scenarios that implement innovations
- breed adjustment
- herd structure
Quantify improvements on farms by comparing with the baseline. Plot improvements vs. the baseline on the above figure. Plot improvements vs. baseline in the emissions intensity figure.
For each slice, and each level of confidence, calculate the percentage of farms that made changes that would be correctly rewarded vs. the ones that did not change (the baseline population). Make a figure showing the Type I and Type II errors resulting from each confidence level at each baseline yield.
Testing interventions
Single-run function
To clearly characterize the impact of interventions, we should compare their impact with a counterfactual, i.e. with the same farms without the intervention. While this is usually impossible in practice, it is easy with a model. We just have to run the simulation again with identical inputs except for those that are affected by the intervention.
We already simulated emissions for a large population of farms, and we can now apply our interventions to this population by modifying certain parameters of the input tables. What we don’t want now, however, is to randomize all values again, since this would add additional noise to the results, which would mask the signal caused by the interventions.
We now first need a function that can apply our greenhouse gas emissions model to a pre-defined set of input values. We can find syntax that can achieve this in the mcSimulation function. Here it is (with slight modifications):
run_model<-function(x_data,model)
{
varnames<-colnames(x_data)[which(!colnames(x_data)=="Scenario")]
model_function_ <- function(x) {
e <- as.list(sapply(X = varnames, FUN = function(i) as.numeric(unlist(x[i]))))
eval(expr = body(model), envir = e)
}
y <- do.call(what = rbind, args = lapply(X = apply(X = x_data,
MARGIN = 1, FUN = model_function_), FUN = unlist))
returnObject <- list(y = data.frame(y), x = data.frame(x_data))
returnObject$call <- match.call()
class(returnObject) <- cbind("mcSimulation", class(returnObject))
return(returnObject)
}The following code can precisely reproduce the results of the previous analysis:
res<-run_model(GHG_simulation_scenarios$x,ghg_emissions)To simulate interventions, we can now feed this function with modified x_data data.frames.
Breed intervention
Introduction of higher-performance breeds may reduce greenhouse gas emissions. We simulate an intervention that replaces animals of low-performance breeds with animals of high-performance breeds.
To retain the farm-specific performance levels obtained from the farm survey, we simulate these impacts in a relative manner. For this, we first need to determine the performance and characteristics of each of the breeds. We’ll base this on the farm survey data (Kenyaherds dataset from above).
cows<-Kenyaherds[which(Kenyaherds$animal_type %in% 4:6),]
ggplot(cows, aes(factor(breed),milk_yield)) +
geom_violin(draw_quantiles=c(0.5),fill="light green") +
xlab("Breed") + ylab("Milk yield (kg per day)") + theme_bw()## Warning: Groups with fewer than two data points have been dropped.
milk_summary<-cbind(aggregate(cows$milk_yield, by=list(cows$breed), FUN=median),table(cows$breed))
milk_summary<-milk_summary[,c(1,4,2)]
colnames(milk_summary)<-c("Breed","n","Median_milk_yield")
milk_summary## Breed n Median_milk_yield
## 1 1 488 6.708090
## 2 2 134 4.999040
## 3 3 16 6.372821
## 4 4 20 5.538156
## 5 5 36 3.035786
## 6 7 1 7.217000
## 7 8 7 3.507483This evaluation shows that the difference in median milk yields is relatively small, but that the milk yield distributions show a lot of variation that is apparently not explained by the breed. It is also apparent that the number of animals these estimates are based on varies considerably. Two breeds make up 88.6 % of the entire population, with others having very few animals, which means that the medians may be somewhat unreliable (especially for the single animal of breed 7).
Let’s define the intervention as follows:
All animals are replaced by high-performance animals of breed 1, which increases milk yields proportionally to the median milk yield differences. This means that milk yields are multiplied by the following factors:
milk_summary[,"Intervention_factor"]<-sapply(milk_summary$Median_milk_yield, function(x) min(x/milk_summary$Median_milk_yield[1],1))
milk_summary ## Breed n Median_milk_yield Intervention_factor
## 1 1 488 6.708090 1.0000000
## 2 2 134 4.999040 0.7452256
## 3 3 16 6.372821 0.9500201
## 4 4 20 5.538156 0.8255936
## 5 5 36 3.035786 0.4525559
## 6 7 1 7.217000 1.0000000
## 7 8 7 3.507483 0.5228735A change of the breed also causes changes in weight, mature weight and weight gain.
live_weights<-aggregate(Kenyaherds$live_weight_LW, by=list(Kenyaherds$animal_type,Kenyaherds$breed), FUN=median)
colnames(live_weights)<-c("animal_type","breed","median_live_weight")
live_weights[,"Intervention_factor"]<-NA
for(i in 1:nrow(live_weights))
live_weights$Intervention_factor[i]<-live_weights$median_live_weight[i]/
live_weights$median_live_weight[which(live_weights$animal_type==live_weights$animal_type[i] & live_weights$breed==1)]
# mature weight should probably be changed, and the code below does this. However, the mature weight is currently the same for all the breeds, so this accomplished nothing.
mature_weights<-aggregate(Kenyaherds$mature_weight_MW, by=list(Kenyaherds$animal_type,Kenyaherds$breed), FUN=median)
colnames(mature_weights)<-c("animal_type","breed","median_mature_weight")
mature_weights[,"Intervention_factor"]<-NA
for(i in 1:nrow(mature_weights))
mature_weights$Intervention_factor[i]<-mature_weights$median_mature_weight[i]/
mature_weights$median_mature_weight[which(mature_weights$animal_type==mature_weights$animal_type[i] & mature_weights$breed==1)]
# same for the weight gain - doesn't depend on the breed at the moment, so all the factors are 1
weight_gains<-aggregate(Kenyaherds$weight_gain_WG, by=list(Kenyaherds$animal_type,Kenyaherds$breed), FUN=median)
colnames(weight_gains)<-c("animal_type","breed","median_weight_gain")
weight_gains[,"Intervention_factor"]<-NA
for(i in 1:nrow(weight_gains))
weight_gains$Intervention_factor[i]<-weight_gains$median_weight_gain[i]/
weight_gains$median_weight_gain[which(weight_gains$animal_type==weight_gains$animal_type[i] & weight_gains$breed==1)]
weight_gains$Intervention_factor[which(weight_gains$median_weight_gain==0)]<-1Now we can modify the input table accordingly.
We go through all the farms and determine the ‘improvement’ coefficients (based on herd structure).
farms<-unique(Kenyaherds$QQNO)
breed_intervention<-data.frame(Farm=farms, milk_yield_4=1, milk_yield_5=1, milk_yield_6=1)
Kenyaherds_interventions<-Kenyaherds[,c("QQNO",
"system",
"animal_type",
"breed",
"live_weight_LW",
"mature_weight_MW",
"weight_gain_WG",
"milk_yield",
"dig_energy",
"CP.")]
Kenyaherds_interventions[,"milk_change"]<-0
cows<-which(Kenyaherds_interventions$animal_type %in% 4:6)
Kenyaherds_interventions$milk_change[cows]<- sapply(cows,function(x) 1/milk_summary$Intervention_factor[milk_summary$Breed==Kenyaherds_interventions$breed[x]] )
milk_change<-aggregate(Kenyaherds_interventions$milk_change[cows], by=list(Kenyaherds_interventions$QQNO[cows]), FUN=mean)
breed_intervention$milk_yield_4[which(breed_intervention$Farm %in% milk_change[,1])]<-milk_change[,2]
breed_intervention$milk_yield_5<-breed_intervention$milk_yield_4
breed_intervention$milk_yield_6<-breed_intervention$milk_yield_4
for(i in 1:nrow(Kenyaherds_interventions))
{Kenyaherds_interventions[i,"live_weight_change"]<-
1/live_weights$Intervention_factor[which(live_weights$animal_type==Kenyaherds_interventions$animal_type[i]&
live_weights$breed==Kenyaherds_interventions$breed[i])]
# Kenyaherds_interventions[i,"mature_weight_change"]<-
# 1/mature_weights$Intervention_factor[which(mature_weights$animal_type==Kenyaherds_interventions$animal_type[i]&
# mature_weights$breed==Kenyaherds_interventions$breed[i])]
Kenyaherds_interventions[i,"weight_gain_change"]<-
1/weight_gains$Intervention_factor[which(weight_gains$animal_type==Kenyaherds_interventions$animal_type[i]&
weight_gains$breed==Kenyaherds_interventions$breed[i])]
}
live_weight_agg<-aggregate(Kenyaherds_interventions$live_weight_change,
by=list(Kenyaherds_interventions$QQNO,Kenyaherds_interventions$animal_type), FUN=mean)
#mature_weight_agg<-aggregate(Kenyaherds_interventions$mature_weight_change,
#by=list(Kenyaherds_interventions$QQNO,Kenyaherds_interventions$animal_type), FUN=mean)
weight_gain_agg<-aggregate(Kenyaherds_interventions$weight_gain_change,
by=list(Kenyaherds_interventions$QQNO,Kenyaherds_interventions$animal_type), FUN=mean)
colnames(live_weight_agg)<-c("farm","animal_type","conv_factor")
#colnames(mature_weight_agg)<-c("farm","animal_type","conv_factor")
colnames(weight_gain_agg)<-c("farm","animal_type","conv_factor")
for(i in 1:nrow(live_weight_agg))
breed_intervention[which(breed_intervention$Farm==live_weight_agg$farm[i]),
paste0("live_weight_LW_",live_weight_agg$animal_type[i])]<-
live_weight_agg$conv_factor[i]
#for(i in 1:nrow(mature_weight_agg))
# breed_intervention[which(breed_intervention$Farm==mature_weight_agg$farm[i]),
# paste0("mature_weight_LW_",mature_weight_agg$animal_type[i])]<-
# mature_weight_agg$conv_factor[i]
for(i in 1:nrow(weight_gain_agg))
breed_intervention[which(breed_intervention$Farm==weight_gain_agg$farm[i]),
paste0("weight_gain_WG_",weight_gain_agg$animal_type[i])]<-
weight_gain_agg$conv_factor[i]
breed_intervention[is.na(breed_intervention)]<-1Go through the x file, for each row take a random draw based on adoption rate (does farm adopt or not), then (possibly) apply factors. Actually, in this case we want to simulate full adoption, so we’ll set this to 1.
x_breed_intervention <- GHG_simulation_scenarios$x
adoption_rate <- 1
for (i in 1:nrow(x_breed_intervention))
{
if (rbinom(1, 1, adoption_rate)) # if farmers adopt, make changes in x file
x_breed_intervention[i, colnames(breed_intervention)[2:ncol(breed_intervention)]] <-
x_breed_intervention[i, colnames(breed_intervention)[2:ncol(breed_intervention)]] *
breed_intervention[which(paste0("Farm_", breed_intervention$Farm) ==
x_breed_intervention$Scenario[i]),
2:ncol(breed_intervention)]
}
# run model for modified x file
breed_res <- run_model(x_breed_intervention, ghg_emissions)Compare outcomes of breed intervention with baseline
breed_res$y[,"CO2eq_per_milk"]<-breed_res$y$pc_on_farm/breed_res$y$pc_milk_yield
breed_res$y<-breed_res$y[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]
breed_res$x<-breed_res$x[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]
impact_of_breed<-data.frame(baseline=dairy_households$y$pc_milk_yield,breed_milk=breed_res$y$pc_milk_yield,CO2_intensity_baseline= dairy_households$y$CO2eq_per_milk,CO2_intensity_breed=breed_res$y$CO2eq_per_milk)
impact_of_breed[,"milk_change"]<-round(impact_of_breed$breed_milk-impact_of_breed$baseline,2)
impact_of_breed[,"emissions_intensity_change"]<-round(impact_of_breed$CO2_intensity_breed-impact_of_breed$CO2_intensity_baseline ,3)
ggplot(data=impact_of_breed,aes(x=milk_change,y=emissions_intensity_change)) + geom_point() + theme_bw() + ylab("Change in emissions intensity (kg CO2-eq kg-1)") + xlab("Change in per capita milk yield (kg)")
ggsave("Fig_XX_breed_reductions.png")## Saving 7 x 5 in imageggplot(data = decreased_emissions_probability, aes(Base_yield, yield_increase)) + geom_contour_filled(aes(z=probability)) +
stat_contour(breaks = c(0.9),aes(z=probability),col="black") +
geom_point(data=impact_of_breed[which(impact_of_breed$milk_change>0),],aes(x=baseline,y=milk_change)) +
ggtitle("Probability that emission intensity reductions occurred") +
xlab("Baseline milk yield per head and year (kg)") +
ylab("Milk yield increase per head and year (kg)") ## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf
ggsave("Fig_XX1_breed_probability_of_intensity_reduction.png")## Saving 7 x 5 in image## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -InfThis plot shows a lot of linear patterns that result from the application of fixed factors to convert the milk yield between cows of different breeds. Not so pretty, but maybe ok.
Herd structure
Retiring unproductive bulls
To simulate this scenario, we simply have to remove all animals of type 1 from the herds.
x_nobull_intervention <- GHG_simulation_scenarios$x
nobull_adoption_rate <- 1.0
for (i in 1:nrow(x_nobull_intervention))
{
if (rbinom(1, 1, nobull_adoption_rate)) # if farmers adopt, make changes in x file
x_nobull_intervention$N_animal_type_1[i]<-0
}
# run model for modified x file
nobull_res <- run_model(x_nobull_intervention, ghg_emissions)Evaluate impact of “Retiring unproductive bulls” intervention
We’re assuming that all bulls >36 months old have been removed.
uncertainty::varscatter(in_var = nobull_res$y$pc_milk_yield,
out_var = nobull_res$y$pc_on_farm,
xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1}))) ## [1] "Scatter plot of influencing variable (x) and outcome variable (y) with associated and estimated densities."uncertainty::varscatter(in_var = GHG_simulation_scenarios$y$pc_milk_yield,
out_var = GHG_simulation_scenarios$y$pc_on_farm,
xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1})))## [1] "Scatter plot of influencing variable (x) and outcome variable (y) with associated and estimated densities."uncertainty::varkernel(in_var = nobull_res$y$pc_milk_yield,
out_var = nobull_res$y$pc_on_farm,
xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1})))## [1] "Density surface plot of an expected influential variable (x) and an outcome variable of interest (y)."uncertainty::varkernel(in_var = GHG_simulation_scenarios$y$pc_milk_yield,
out_var = GHG_simulation_scenarios$y$pc_on_farm,
xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1})))## [1] "Density surface plot of an expected influential variable (x) and an outcome variable of interest (y)."
No bulls intervention (left) vs. baseline (right)
Reducing unproductive replacement males
The idea is to raise the off-take rate of replacement males by 23% (not sure why this number). We’ll implement this by eliminating replacement males on all adopting farms, with a 23% adoption rate. The result should be quite similar.
x_lessReplace_intervention <- GHG_simulation_scenarios$x
lessReplace_adoption_rate <- 0.23
for (i in 1:nrow(x_lessReplace_intervention))
{
if (rbinom(1, 1, lessReplace_adoption_rate)) # if farmers adopt, make changes in x file
x_lessReplace_intervention$N_animal_type_3[i]<-0
}
# run model for modified x file
lessReplace_res <- run_model(x_lessReplace_intervention, ghg_emissions)Compare outcomes of herd intervention with baseline
No bulls scenario
nobull_res$y[,"CO2eq_per_milk"]<-nobull_res$y$pc_on_farm/nobull_res$y$pc_milk_yield
nobull_res$y<-nobull_res$y[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]
nobull_res$x<-nobull_res$x[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]
impact_of_nobull<-data.frame(baseline=dairy_households$y$pc_milk_yield,breed_milk=nobull_res$y$pc_milk_yield,CO2_intensity_baseline= dairy_households$y$CO2eq_per_milk,CO2_intensity_breed=nobull_res$y$CO2eq_per_milk)
impact_of_nobull[,"milk_change"]<-round(impact_of_nobull$breed_milk-impact_of_nobull$baseline,2)
impact_of_nobull[,"emissions_intensity_change"]<-round(impact_of_nobull$CO2_intensity_breed-impact_of_nobull$CO2_intensity_baseline ,3)
ggplot(data=impact_of_nobull,aes(x=milk_change,y=emissions_intensity_change)) + geom_point() + theme_bw() + ylab("Change in emissions intensity (kg CO2-eq kg-1)") + xlab("Change in per capita milk yield (kg)")
ggplot(data = decreased_emissions_probability, aes(Base_yield, yield_increase)) + geom_contour_filled(aes(z=probability)) +
stat_contour(breaks = c(0.9),aes(z=probability),col="black") +
geom_point(data=impact_of_nobull[which(impact_of_nobull$milk_change>0),],aes(x=baseline,y=milk_change)) +
ggtitle("Probability that emission intensity reductions occurred") +
xlab("Baseline milk yield per head and year (kg)") +
ylab("Milk yield increase per head and year (kg)") ## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf
ggsave("Fig_XX1_breed_probability_of_intensity_reduction.png")## Saving 7 x 5 in image## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Infggsave("Fig_XX_impact_of_nobull_reductions.png")## Saving 7 x 5 in image## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -InfFewer replacement males scenario
lessReplace_res$y[,"CO2eq_per_milk"]<-lessReplace_res$y$pc_on_farm/lessReplace_res$y$pc_milk_yield
lessReplace_res$y<-lessReplace_res$y[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]
lessReplace_res$x<-lessReplace_res$x[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]
impact_of_lessReplace<-data.frame(baseline=dairy_households$y$pc_milk_yield,breed_milk=lessReplace_res$y$pc_milk_yield,CO2_intensity_baseline= dairy_households$y$CO2eq_per_milk,CO2_intensity_breed=lessReplace_res$y$CO2eq_per_milk)
impact_of_lessReplace[,"milk_change"]<-round(impact_of_lessReplace$breed_milk-impact_of_lessReplace$baseline,2)
impact_of_lessReplace[,"emissions_intensity_change"]<-round(impact_of_lessReplace$CO2_intensity_breed-impact_of_lessReplace$CO2_intensity_baseline ,3)
ggplot(data=impact_of_lessReplace,aes(x=milk_change,y=emissions_intensity_change)) + geom_point() + theme_bw() + ylab("Change in emissions intensity (kg CO2-eq kg-1)") + xlab("Change in per capita milk yield (kg)")
ggplot(data = decreased_emissions_probability, aes(Base_yield, yield_increase)) + geom_contour_filled(aes(z=probability)) +
stat_contour(breaks = c(0.9),aes(z=probability),col="black") +
geom_point(data=impact_of_lessReplace[which(impact_of_lessReplace$milk_change>0),],aes(x=baseline,y=milk_change)) +
ggtitle("Probability that emission intensity reductions occurred") +
xlab("Baseline milk yield per head and year (kg)") +
ylab("Milk yield increase per head and year (kg)") ## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf
ggsave("Fig_XX1_breed_probability_of_intensity_reduction.png")## Saving 7 x 5 in image## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Infggsave("Fig_XX_lessReplace_reductions.png")## Saving 7 x 5 in image## Warning: stat_contour(): Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf